function U = poisson_2d(left, right, bottom, top, Nx, Ny, basis_type, N_gauss_int2d, function_cor, function_f, function_g, boundary_type)

[node, elem, bdary] = processing_node_elem_bdary_2d(left, right, bottom, top, Nx, Ny, basis_type, boundary_type);

dN = reference_basis_function_2d(basis_type);

A1 = assemble_matrix_A_2d(node, elem, function_cor, dN, basis_type, "dx", "dx", N_gauss_int2d);

A2 = assemble_matrix_A_2d(node, elem, function_cor, dN, basis_type, "dy", "dy", N_gauss_int2d);

A = A1 + A2;

b = assemble_vector_b_2d(node, elem, function_f, dN, basis_type, "x", N_gauss_int2d);

[A, b] = processing_boundary_condition_2d(A, b, node, bdary, function_g);

U = A\b;

end